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Abstract - We study a system formed by soft colloidal spheres attracting each other via a 
square-well potential, using extensive Monte Carlo simulations of various nature. The softness 
is implemented through a reduction of the infinite part of the repulsive potential to a finite one. 
For sufficiently low values of the penetrability parameter we find the system to be Ruelle stable 
with square- well like behavior. For high values of the penetrability the system is thermodynami- 
cally unstable and collapses into an isolated blob formed by a few clusters each containing many 
overlapping particles. For intermediate values of the penetrability the system has a rich phase 
diagram with a partial lack of thermodynamic consistency. 



, Pair effective interactions in soft condensed matter 
physics can be of various nature and one can often find real 
systems whose interaction is bounded at small separations 
as, for instance, in the case of star and chain polymers 
[J]. In this case, paradigmatic models, such as square-well 
('SW) fluids, that have been rather successful in predicting 
lihermo-physical properties of simple liquids, are no longer 
useful. Instead, different minimal models accounting for 
the boundness of the potential have to be considered, the 
Gaussian core model [5] and the penetrable-sphere (PS) 
ipodel [5H5] being well studied examples. More recently, 
the penetrable square- well (PSW) fluid has been added to 
this category |BHH] with the aim of including the existence 
of attractive effective potentials. The PSW model is ob- 
tained from the SW potential by reducing to a finite value 
the infinite repulsion at short range. 



4>psw{r) = 




(1) 



where and are two positive energies accounting for 
the repulsive and attractive parts of the potential, respec- 
tively, A is the width of the attractive square well, and 



a is the width of the repulsive barrier. For — oo one 
recovers the SW model, while for A = or = one 
recovers the PS model. 

For finite e^, ea/^r is a measure of the penetrability of 
the barrier and we shall refer to Ca/cr as the penetrability 
ratio. PSW pair potentials can be obtained as effective po- 
tentials for instance in polymer mixtures pUllTTj . While in 
the majority of the cases the well depth ea is much smaller 
than the repulsive barrier (low penetrability limit) this 
mesoscopic objects are highly sensitive to external condi- 
tions (e.g. quality of the solvent) and may thus in principle 
exhibit higher values of the penetrability ratio ea/sr- 

It is well known that three-dimensional SW fluids ex- 
hibit a fluid-fluid phase transition for any width of the 
attractive square well |12H16| . the liquid phase becoming 
metastable against the formation of the solid for a suffi- 
ciently narrow well [TS]. It is also well established that 
in the PS fluid (that lacking an attractive component in 
the pair potential cannot have a fluid-fluid transition) an 
increase of the density leads to the formation of clusters 
of overlapping particles arranged in an ordered crystalline 
phase [linHin]. 

While the novel features appearing even in the one- 
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dimensional case have been studied in some details [SHS], 
no analysis regarding the influence of penetrability and 
attractiveness on the phase behavior of PSW fluids have 
been reported, so far, in three dimensions. The present 
paper aims to fill this vacancy. 

A system is defined to be Ruelle stable when the to- 
tal potential energy Un for a system of N particles sat- 
isfies the condition Un > —NB where i? is a finite pos- 
itive constant [20l[2T]. In Ref. [8] we proved that in the 
one-dimensional case the PSW model is RucUc stable if 
ea/cr < 1/2(^-1- 1), where £ is the integer part of A/a. 
This result can be extended to any dimensionality d by 
the following arguments. The configuration which mini- 
mizes the energy of N particles interacting via the PSW 
potential is realized when M closed packed clusters, each 
consisting of s = N/M particles collapsed into one point, 
are distributed such that the distance between centers of 
two neighbor clusters is a. In such a configuration, all 
the particles of the same cluster interact repulsively, so 
the repulsive contribution to the total potential energy 
is erMs{s — l)/2. In addition, the particles of a given 
cluster interact attractively with all the particles of those 
/a clusters within a distance smaller that a + A. In the 
two-dimensional case, /a = 6 and 12 if A/a < V^- 1 and 
•^3— 1 < A/ct < 1, respectively. For d = 3, the case we are 
interested in, one has /a = 12, 18, and 42 HA/a < a/2-1, 
\/2 - 1 < A/a < \/3 - 1, and \/3 - 1 < A/a < 1, respec- 
tively. The attractive contribution to the total potential 
energy is thus -ea(M/2) [/a - bA{M)] s^, where 6a(M) 
accounts for a reduction of the actual number of clus- 
ters interacting attractively, due to boundary effects. This 
quantity has the properties 6a (Af) < /a, 6a (1) = /a, and 
limj\f_j.oo b\[M) = 0. For instance, in the two-dimensional 
case with A/(7 < VS-lone has 6a(M) = 2{Ay/M-l)/M. 
Therefore, the total potential energy is 



Un{M) 
Ner 



1 N 

— + F(M), 

2 2M ^ ^' 



(2) 



where F(M) = (ea/e,)6A(M) + (1 - fA^a/er)- If ea/e^ < 
1//a, F{M) is positive definite, so Un{M)/N has a lower 
bound and the system is stable in the thermodynamic 
limit. On the other hand. If ea/er > 1//a, one has 

F(l) = 1 but limA/^oo = -(/AEa/er " 1) < 0. In 

that case, there must exist a certain finite value M = Mq 
such that F{M) < for M > Mq. As a consequence, 
in those configurations with M > Mq, Un{M)/N has no 
lower bound in the limit — > oo and thus the system may 
be unstable. 

We have performed an extensive analysis of the vapor- 
liquid phase transition of the system using Gibbs Ensem- 
ble Monte Carlo (GEMC) simulations [HH^^, starting 
from the corresponding SW fiuid condition and gradually 
increasing the penetrability ratio Ca/cr until the transition 
disappears. A total number of = 512 particles with 2N 
particle random displacement, N/10 volume changes, and 
N particle swap moves between the gas and the liquid 
box, on average per cycle, were considered. We find that 



for any given width A/ a < 1 of the well, there is a limit 
value of the penetrability ratio ea/sr above which no fluid- 
fluid phase transition is observed. This is depicted in Fig. 
[T] where it can be observed that (for A/ a < 1) this line 
lies outside the Ruelle stable region Ca/er < 1//a- 
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Fig. 1: Plot of the penetrability ratio e^/er as a function of 
A/cr. The displayed line separates the parameter region where 
the PSW model admits a fluid-fluid phase transition from that 
where it does not. The highlighted region (ea/tr < 1/12 for 
A/a < a/2 - 1, ea/e,. < 1/18 for a/2 - 1 < A/a < a/3 - 1, 
and Ea/er < 1/42 for a/3 — 1 < A/a < 1) shows where the 
model is expected to be thermodynamically stable in the sense 
of Ruelle for any thermodynamic state. The SW model falls 
on the horizontal axis (ea/tr — >■ 0) and its fluid-fluid transition 
is expected to be metastable against the freezing transition 
for A/a < 0.25 [inj. The circles are the points chosen for 
the calculation of the coexistence lines (see Figs. [2]and[5|. 
The crosses are the points chosen for the determination of the 
boundary between extensive and non-extensive phases (see Fig. 

It is instructive to analyze the detailed form of the 
coexistence curves below (but close to) the limit line of 
Fig. [TJ This is presented in Fig. [5J We have explicitly 
checked that our code reproduces completely the results 
of Vega et al. [H] for the SW model. Following stan- 
dard procedures [12] we fitted the GEMC points, near 
the critical point using the law of rectilinear diameters 
{pi + pg)/2 ^ pc + A{Tc - T) where pi [pg) is the density 
of the liquid (gas) phase, pc is the critical density, and 
Tc is the critical temperature. Furthermore, the tempera- 
ture dependence of the density difference of the coexisting 
phases is fitted to the scaling law pi — pg = B{Tc — T)^ 
where (3 = 0.32 is the critical exponent for the three- 
dimensional Ising model. The amplitudes A and B where 
determined from the fit. In the state points above the 
limit line of Fig. [T]we have considered temperatures below 
the critical temperature of the corresponding SW system. 
The disappearance of the fluid-fluid transition is signaled 
by the evolution towards an empty gas box and a clustered 
phase in the liquid box. 

As discussed, the PSW fiuid is thermodynamically Ru- 
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such that the system is metastable if T > Tinst and unsta- 
ble if T < Tinst- We determined the metastable/unstable 
crossover by performing NVT Monte Carlo simulations 
with N = 512 particles initially uniformly distributed 
within the simulation box. Figure [3] reports the results 
in the reduced temperature-density plane for A/cr = 0.5 
and for two selected penetrability ratios Sa/^r = 1/7 and 
1/4, the first one lying exactly just above the limit line of 
Fig. [T] while the second deep in the non-transition region. 
We worked with constant size moves (instead of fixing the 
acceptance ratios) during the simulation run, choosing the 
move size of 0.15 in units of the the simulation box side. A 
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Fig. 2: Fluid-fluid coexistence line. The solid circle represents 
the critical point (pc, Tc). In the top panel, for A/a = 0.25 and 
Ea/fr ~ 1/6, one has pcO-^ ~ 0.307 and ksTc/ea = 0.762; in the 
middle panel, below the limit penetrability, for A/a = 0.5 and 
Sa/sr = 1/8, one has pca^ — 0.307 and ksTc/ea = 1.241; and 
in the bottom panel for A/a — 1.0 and ta/er = 1/11, one has 
Pca^ = 0.292 and ksT^/ea = 2.803. The lines are the result of 
the fit with the law of rectilinear diameters. The SW results 
are the ones of Vega et al. [T5] 

elle stable when ea/^r < 1//a for all values of the ther- 
modynamic parameters. For Ca/cr > 1//a the system is 
either extensive or non-extensive depending on temper- 
ature and density. For a given density, one could then 
expect that there exists a certain temperature Tinst (p), 



Fig. 3: Regions of the phase diagram where the PSW fluid, 
with A/a = 0.5 and two different values of ta/^r, is found to 
be extensive or non-extensive (here we used A'^ = 512 particles). 
The left shows representative snapshots in the two regions. 
The instability line corresponding to the higher penetrability 
case (dashed line) lies above the one corresponding to lower 
penetrability (continuous line). 

crucial point in the above numerical analysis is the iden- 
tification of the onset of the instability. Clearly the phys- 
ical origin of this instability stems from the fact that the 
attractive contribution increases unbounded compared to 
the repulsive one and particles tend to lump up into clus- 
ters of multiply overlapping particles ("blob"). Hence the 
energy can no longer scale linearly with the total number 
of particles N and the thermodynamic limit is not well 
defined (non-extensivity) . We define a cluster in the fol- 
lowing way. Two particles belong to the same cluster if 
there is a path connecting them, where we are allowed to 
move on a path going from one particle to another if the 
centers of the two particles are at a distance less than a. 

The state points belonging to the unstable region are 
characterized by a sudden drop of the internal energy and 
of the acceptance ratios at some points in the system evo- 
lution during the MC simulation. Representative snap- 
shots show that a blob structure has nucleated around a 
certain point and occupies only a part of the simulation 
box with a few clusters. The number of clusters decreases 
as one moves away from the boundary line found in Fig. 
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[3] towards lower temperatures. Upon increasing ta/^r the 
number of clusters decreases and the number of particles 
per cluster increases. We assume a state to be metastable 
if the energy does not undergo the transition after 10^ 
single particle moves. The fluid-fluid transition above the 
limit penetrability line of Fig. [T] is not possible because 
the non-extensive phase shows up before the critical point 
is reached. 

The boundary line of Fig. [3] is robust with respect to 
the size of the system, provided that a sufficiently large 
size {N > 512) is chosen. When the number of particles 
in the simulation goes below the number of clusters which 
would form in the non-extensive phase the system seems to 
remain extensive. For instance, with N = 1024 particles 
we obtained under the ea/^r ~ 1/7, A/ct = 0.5 conditions 
a threshold temperature kBT/ea ~ 1-15 for pa^ ~ 0.4 and 
2.25 for pa^ = 0.8, which are close to the values obtained 
with A^ = 512 particles. 

There is an apparent hysteresis in forming and melting 
the non-extensive phase. For example when Sa/^r = 1/7, 
A/a = 0.5, and pa^ = 1.0 the non-extensive phase starts 
forming when cooling down to ksT/ea = 2.75. Upon in- 
creasing the temperature again, we observed a melting 
transition at significantly higher temperatures {kBT/ea > 
4). We also found the hysteresis to be size dependent; in 
the same state for pa^ = 0.6 the melting temperatures are 
ksT/ea w 2.5 for A^ = 256, ksT/ea ~ 4.5 for A = 512, 
and ksT/ea « 6.5 for A^ = 1024. This suggests that the 
extensive phase in Fig. |2] is actually metastable with re- 
spect to the non-extensive phase in the thermodynamic 
limit. However the metastable phase can be stabilized by 
taking the size of the system finite. In addition we cannot 
exclude, a priori, the possibility of a true extensive stable 
phase as it is not prevented by the Ruelle criterion. We 
note that the size dependence of the hysteresis in the melt- 
ing could be attributed to the fact that the blob occupies 
only part of the simulation box and therefore a surface 
term has a rather high impact on the melting tempera- 
ture. 

A convenient way to characterize the structure of the 
fluid is to consider the radial distribution function g{r). 
This is depicted in Fig. |4] for the cases A/cr = 0.5, 
ksT/ta = 1-20, and pa^ = 0.7 at ea/^r = 1/8 and 
Ea/fir = 1/7. The latter case is in the non-extensive re- 
gion, according to Fig. [31 We can clearly see that there 
is a dramatic change in the structural properties of the 
PSW hquid. In the non-extensive case, ea/e^ = 1/7, the 
radial distribution function grows a huge peak at r = 
and decays to zero after the first few peaks, which sug- 
gests clustering and confinement of the system. 

In order to study the solid phase of the PSW model 
below the limit penetrability we employed isothermal- 
isobaric (NPT) MC simulations. A typical run would con- 
sist of 10* steps (particle moves or volume moves) with 
an equilibration time of 10'' steps. We used 108 particles 
and adjusted the particle moves to have acceptance ratios 
of w 0.5 and volume changes to have acceptance ratios 




Fig. 4: Radial distribution function for the PSW model at 
A/o- = 0.5, ksT/ea = 1.20, and pa^ = 0.7 for two different 
values of the penetrability ratio ea/tr- 

of w 0.1. Here we only consider the case of PSW with 
A/cr = 0.5 and ta/^r = 1/8 and 1/15. 

For the SW system with a width A/cr = 0.5 the 
critical point is known to be at ksTc/ta ~ 1-23 and 
Pccr^ — 0.309, its triple point being at ksTt/ea — 0.508, 
Pta^/ea ^ 0.00003, pia^ = 0.835, and p^cr^ = 1.28 [H]. 
No solid stable phase was found in Ref. [12] for temper- 
atures above the triple point, meaning that the melting 
curve in the pressure-temperature phase diagram is al- 
most vertical. On the other hand, the phase diagram of 
the PSW fluid with the same well width and a value of 
Ca/cr = 1/8, just bclow the limit line of Fig. [TJ shows 
that the melting curve has a smooth positive slope in the 
pressure-temperature phase diagram. In order to establish 
this, we used NPT simulations to follow the fcsT/ea = 1 
isotherm. From Fig. [5]we can clearly sec the jumps in den- 
sity corresponding to the gas-liquid and to the liquid-solid 
coexistence regions. The presence of a solid phase can be 
checked by computing the Qe order parameter |27] , calcu- 
lated for the center of mass of individual clusters, that in 
the present case turns out to be Qq « 0.35. The crystal 
structure is triclinic with a unit cell with a = 6 = c = ct 
and a = (3 — tt/3 and 7 = cos~^(l/4). There are possibly 
other solid-solid coexistence regions at higher pressures. 
Moreover the relaxation time of the MC run in the solid 
region is an order of magnitude higher than the one in the 
liquid region. 

We also run at the same temperature a set of simulations 
for the PSW fluid with Ca/cr = 1/15. The results (see Fig. 

showed no indication of a stable solid, in agreement with 
the fact that at this very low value of the penetrability 
ratio the system is SW-like. 

A peculiarity of the PSW in the region below the limit 
penetrability of Fig. [U but not in the Ruelle stability re- 
gion, is a violation of the Clausius-Clapeyron equation [5S] 
along the liquid-solid coexistence curve, what represents a 
partial lack of thermodynamic consistency 
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Fig. 5: Isotherm ksT/ea = 1 for the PSW system with A/<t = 
0.5 and ta/^r ~ 1/8 and ea/tr ~ 1/15, as obtained from NPT 
MC simulation with A'^ = 108 particles. The pressure axis is in 
logarithmic scale. The right shows snapshots of the centers of 
mass of the clusters in the solid. 

In the intermediate penetrability case (i.e., above Ru- 
elle's threshold but below the limit penetrability), the ob- 
served crystal structure is made of clusters of overlapping 
particles (rarely more than two) located at the sites of a 
regular crystal lattice. It is precisely this additional de- 
gree of penetrability, not present in the SW system, that 
allows for the coexistence of the liquid and the solid at 
not excessively large pressures. In this respect qualitative 
arguments along the lines suggested in Ref. [T^ could be 
useful. 

In conclusion, we have studied the phase diagram of 
the three-dimensional PSW system. This model is Ru- 
elle stable for ea/cr < 1//a- For Ca/cr > 1//a is either 
metastable or unstable (non-extensive), depending on the 
values of temperature and density, as shown in Fig. [31 The 
instability is indicated by the collapse of the system in a 
confined blob made up of a few clusters of several overlap- 
ping particles. Moreover, the gas-liquid phase transition 
disappears, as shown in Fig. [TJ 

For the metastable fluid near the limit penetrability line 
of Fig. [T] we determined the phase diagram comparing it 
with the corresponding SW case. We determined how the 
gas-liquid coexistence curves are modified by the presence 
of penetrability (see Fig. [5]) and discussed the main fea- 
tures of the phase diagram, including the solid phase, for 
A/a = 0.5. 

For the liquid-solid coexistence curves we generally 
found that the solid density increases with respect to the 
corresponding SW case, as expected, due to the forma- 
tion of clusters of overlapping particles in the crystal sites. 
For A/cr = 0.5 the PSW model with a sufficient pene- 
trability to have a metastable system, but not a Ruelle 
stable one, has a melting curve with a positive slope in 
the pressure-temperature phase diagram with a violation 
of the Clausius-Clapeyron thermodynamic equation, thus 
confirming the metastable character of the phases. For 
sufficient low penetrability the system is in the Ruelle sta- 



ble region, and behaves as the corresponding SW model. 

In summary, by experimentally tuning the repulsive bar- 
rier relative to the well depth one could observe (a) stable 
phases resembling those of a normal fluid, (b) metastable 
phases with fluid-fluid and fluid-solid coexistence, or (c) 
the collapse of the system to a small region. 
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